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METHODOLOGY  FOR  STOCHASTIC  MODELING 


1.  INTRODUCTION 

The  requirement  to  develop  stochastic  mathematical  models  arises  across 
the  whole  range  of  engineering  and  applied  research  where  observations  are  made 
of  a  physical  process,  corrupted  by  noise,  and  It  Is  desired  to  determine  the 
underlying  nature,  either  In  time  or  frequency,  of  the  observed  phenomenon. 

This  paper  presents  some  of  the  current  approaches  In  stochastic  modeling  which 
have  been  found  to  be  useful  In  this  area. 

Consider  a  system  In  which  an  Input  at  produces  an  output  yt  as  shown  in 
Figure  1. 


Input  or 
Excitation 


at 


System 

yt 

Output 


Figure  1 

In  general  we  entertain  the  idea  of  expressing  the  output  yt  in  terms  of 
its  past  history  yt-l,  yt-2>  •••>  Yt-p  an(1  the  current  and  past  history  of  the 
excitation  at,  at-l»  at-2  •••  at-q  as  shown 
q  P 

n  =  l  9kat-k  -  1  4>k  yt  -k  m 

k=0  k=l  m 

with  6q  =1 

where  the  coefficients  9j<  and  $)<  are  to  be  estimated  as  well  as  the  parameters 
p  and  q.  Equation  (1)  represents  an  autoregressive  moving  average  (ARMA)  model 
of  order  (p,q).  It  is  advantageous  to  assume  that  the  input  or  excitation  at 
is  a  zero  mean  gaussian  white  noise.  For  the  case  where  all  <^=0  we  have 
what  is  called  a  moving  average  (MA)  model. 


(2) 


T* 


q 

n  ■  I  9k  «t-k 

k=Q 

and  for  6^  =  0  ( k >1)  we  have  an  autoregressive  (AR)  model 
P 

yt  "  *  I  H  yt-k  +  at  (for  e0  =1) 

k=l 

2.  BOX-JENKINS  TECHNIQUE 

The  Box-Jenkins  [1]  approach  for  time  series  analysis  of  observed  data  is 
widely  used  and  is  based  on  the  statistical  properties  of  the  data,  namely,  the 
autocorrelation  function  (ACF)  and  the  partial  autocorrelation  function  (PAC)  for 
each  of  the  AR  and  MA  models. 

a.  Autoregressive  (AR)  model  of  order  p:  An  autoregressive  model  of 
order  p,  AR(p),  is  characterized  as  follows: 

(1)  The  autocorrelation  function  (ACF)  is  either  a  decreasing  exponential 
or  a  damped  sine  wave. 

(2)  The  partial  autocorrelation  function  (PCF)  is  non  zero  for  lags  less  or 
equal  to  p  and  zero  for  lags  greater  than  p,  i.e.,  the  PCF  cutoff  after  lag  p. 

When  we  multiply  equation  (3)  by  yt_k  (k>0)  and  take  expectation  we  obtain 
the  Yule-Walker  equations 


p|c  -  4>1  Pk-l  •  $2  p k-2  -  •••  -  <f>p  Pk-p 

(k  =  1,  2,  ....  p) 


(4) 


where  Pk  is  the  autocorrelation  function  of  the  process  yt  of  lag  k.  The  Yule- 
Walker  equations  can  be  solved  for  $j  when  we  replace  the  theoretical  auto¬ 
correlation  pk  by  the  estimated  autocorrelation  rk,  i.e.,  where 


-1 

£  =pp  £p 


(5) 
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i.  =Pp  £p 


1 
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P2 
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1  PI  P2  • 
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•  •  Pp-1 
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II 
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■  P  = 
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* 

-  P  - 

_Pp_ 

-Pp-1  Pp-2  Pp-3  • 

•  •  pi  JL 

Since  the  autocorrelation  function  for  an  AR(p)  is  infinite  in  extent,  the 

partial  autocorrelation  $4ji,also  called  "reflection  coefficient,"  can  be  obtained 

from  the  Levinson-Ourbin  Algorithm  [2,  3]  namely 

/  rj  2=1 

2-1 

r  2  "  J.  ^  2-1 , j  r  2-j 


*22  =< 


2-1  - 

1  -  ^  *2-l,j  rJ 


2  =2,  3,  ....  L  (5) 


where 


*2j  =  *2-1 , j  ”  *22  *2-1,  2-j 


( j-1 ,2, . . . ,  2-1.) 


where  rj  is  the  autocorrelation  of  the  sampled  data  of  lag  j  with  the  initial 
values  $21  and  $22  given  by 

n  (l-r2) 

*21  *  - 

1  -  r 


*22  = 


(7) 


The  95  percent  confidence  interval,  or  (+2 a  ),  for  rk  and  $|<ic  are  given  by 


Var  [rk] 


q 

1  +  2  V 
i  =  l 


;  (k  >  q) 
9 


(8) 

( Bartel l's  approximation) 
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which  are  useful  quantities  In  determining  order  of  the  model. 

It  Is  advantageous  to  plot  the  ACF  (r^)  and  PCF  Ukfc)  for  lag  1  up  to  any 
specified  lag  k  with  the  corresponding  95  percent  confidence  Interval  +2  o. 

The  selection  of  the  form  of  the  model  is  determined  from  the  behavior  of  the 
autocorrelation  function  and  the  partial  autocorrelation  function  as  follows: 

a.  If  the  autocorrelation  function  decays  either  exponentially  or  sinuos- 
idally,  and  the  partial  autocorrelation  function  ("reflection  coefficient") 
cuts  off  after  a  lag  p,  then  the  model  is  AR  (p)  i.e.,  no  moving  average  terms 
exist. 

b.  If  the  autocorrelation  function  cuts  off  after  lag  q,  and  the  partial 
autocorrelation  function  decays  either  exponentially  or  sinusoidally,  then  we 
have  a  moving  average  model  MA(q)  i.e.,  no  autogressive  terms  exist. 

c.  Should  there  be  no  cutoff  in  either  the  autocorrelation  function  or 
partial  autocorrelation  function,  and  they  independently  decay  either  exponen¬ 
tially  or  sinusoidally,  then  the  model  is  ARMA  (p,  q)  i.e.,  autoregressive-moving 
average  of  order  p,  q  respecti vely,  so  that  there  exist  both  autoregressive  and 
moving  average  terms. 

The  Box-Jenkins  procedures  for  ARMA  modeling  are  given  as  follows: 

Program  #1 : 

Define  expected  value  and  variar,__'  as  follows: 


(Expected  value) 


(Variance) 

Autocovariance  function: 


1  N 

x  =  —  r  x^ 
N  t=l 

c2  .  r 
sx  '  c0 


2i 


To  estimate  e  of  the  moving  average  parameters  we  use  the  autocovariance 


of  xt,  namely  c^,  and  calculate 

(  P  p  .  *  . 

I  l  $1*0  HO  c  Ij+i-kl  P>(3  (4>00  =  "1) 

i=0  k=0  1  1 

cj  P=0  (j=0,  1 . q) 

Box  [1]  then  applies  to  the  Newton-Raphson  algorithm 


Ti+1  = 


-  h 


where 

T1'] =  _f 1 

To  calculate  the  vectors  x1+*  at  the  (i+1)  iteration  from  its  value  x ’  at 
the  ith  iteration 


1  =  (to,  t i,  ....  xq) 
*  qfj 

fj  Y=0  Ti  Ti+J*  " Cj’ 
r  *  (fO.  fl.  •••  ,  fq) 


T0  T1 . Tq 

TO  Ti  .  .  .  .  Tq 

t  1  x  2  ...  Tq 

• 

•  T0  T1  •  •  Tq-1 

• 

+ 

•  »  • 

•  • 

•  •  • 

• 

•  • 

Tq  o 

o  .  . 

—  — 

—  10. 

Newton-Raphson  algorithm  converges  when  | f -j | <  e,  j=0,  1,  ....  q  for  some 
prescribed  e  with 


and  the  estimate  a\  of  white  noise  variance 

q>0 

q=0 

« 

We  therefore  have  calculated  all  the  parameters  in  the  ARMA  (p,  q)  model. 

3.  OVERDETERMINED  RATIONAL  MODELS 

Cadzow  [4]  has  developed  spectral  estimates  using  an  overdetermined  set  of 
Yule-Walker  equations  for  a  rational  model.  The  procedures  are  sunmarlzed  below: 
1.  For  spectral  models  employing  exact  autocorrelation  lag  information. 

a.  Moving  average  model  (MA),  a  B1 ackman-Tukey  approach  is  used  for 
determining  the  spectrin  Sx  (e1a)  )  namely 

q 

Sx  (eJu)  *  l  w(n)  rx  (n)  e^*1*0  no) 

na-q  -  1 

where  w(n)  is  a  desired  data  window,  a  wide  variety  of  which  can  be  located  in 

[5,  6,  7]. 

rx(n)  is  the  exact  autocorrelation  function  of  lag  n  for  the  observed 
signal  x(t). 

b.  Autoregressive  model  (AR)  -  Having  determined  the  order  p  of  a  purely 
autoregressive  model  from  Box-Oenkins  technique  by  studying  partial  autocorrela¬ 
tions,  we  form  a  (p+l)x(p+l)  AR  autocorrelation  matrix  R  whose  elements  are 

R(i,j)  s  rx  (i  -  j)  U  i  <  p  +  1  (11) 

U  j  <  p  +  1 


{  T§ 


*a*  { 


c0  -  I  *iCi 
i*l 


ap) 


Solve  for  Ra  a  jbQ|2  ej_ 
where  _a  =  (1,  aj,  a  2  .... 


(12) 


u  ■  i,"  i."  v1"  -  ■  i 


and  bp  is  selected  so  that  the  first  component  of  &  is  one. _ _ 

Using  the  value  of  _a  just  calculated,  the  spectrun  sx  («*■  )  is  given  by 


SY  (e>) 


1+a^  e"J“  +  ...  +  ap  e 


-jpw 


(13) 


A  more  efficient  method  for  calculating  recursively  the  components  of  a  and  at 
the  same  time  determining  the  order  of  the  AR  model  is  the  Levinson-Durbin 
algorithm.  The  procedure  for  coefficient  determination  is  as  follows: 

Step  #1 

•1(1)  *  -  rx  (l)/rx(0)  (14) 

l^l2  *  [1  -  |  ap) J2]  rx  (0) 

where  we  define  a^k)  as  the  jth  coefficient  associated  with  the  kth  iteration 
Step  #2.  For  k  3  2,  3,  4  ... 


4k)  =  -  c  rx(k) +  i.  ^  (k-1)  rx  (k‘m)]/  i  b4k_1)i2  (is) 

m=l 

a|k)  =  a|k-1)  +  aj(k)  ak_(k*^  1<  i  <  (k-1) 

l4k)l2  -  C  1  -  I  4k)  I2  ]  Ib^-D  |2 

The  iteration  stops  when  assumes  a  constant  value  or  a  sufficiently  small 
value  after  some  value  of  "k"  which  thus  establishes  the  order  of  the  model. 

c.  Autoregressive  moving  average  (ARMA)  -  In  general  an  ARMA  model  takes 
the  form 

P  q 

x(n)  +  I  a|<  x  (n-k)  ■  J  b|<  e  (n-k) 


(16) 


where  [  e(n)j  is  normalized  white  noise.  For  a  model  to  be  causal  the  Yule- 

Walker  equations  take  on  a  simple  form  for  n>q 

p 

n>  q  +  1 


l  ak  rx  (n-k)  =  0 
k*0 


(17) 


For  an  overdetermined  model,  the  extended  Yule-Walker  equations  for  q+1  £  n  <  q+t 
i.e.,  t  linear  equations  in  p  autoregressive  parameter  unknowns,  is  given  by 


rx  (q+1)  rx  (q) . rx  (q-p+1) 

rx  (q+2)  rx  (q+1)  .  .  .  .  rx  (q-p+2) 

1 

al 

1 

o  o 

1 

a2 

• 

= 

0 

• 

_rx  (q+t)  rx  (q+t-1)  .  .  .  rx  (q-p+t)  _ 

• 

• 

0 

or  simple  Ri  _a  *  0 

where  Rj  is  a  Toeplitz  matrix 

R1  ( 1  * J)  *  rx  (q  +  1  +  i  -  j)  1<  i  <  t 

1<  j  <  p  +  1 

where 

±  -  (1»  al t  a2 »  •••  ap) 

are  the  (p+1)  autoregressive  coefficients  and 

rx  (n)  =  E  {  x  (n  +  m)  7  (m)} 

Using  the  above  developed  structure  the  algorithm  for  ARMA  model  is  as  follows 

1.  Form  the  t  x  (p+1)  ARMA  autocorrel ati on  matrix  Rj.  Define  R^  as  the 
transpose  complex  conjugate  matrix. 

2.  If  the  rank  (R^  Rj)  <  p+1  then  solve 

R1  hi  =  o 

using  singular  value  decomposition  (SVD)  to  determine  "p”.  Discussion  on  SVD 
will  be  presented  later  in  this  paper. 


3.  If  the  rank  (Rj  Rj)  *  p+1  then  we  solve 

*1  h  «  ■  «£i  (19) 

where  ej  a  (1,  0,  ...0)'  and  a  Is  a  normalizing  factor  so  that  the  first 

component  of  £  Is  one  or  we  can  find  the  minimum  eigenvalue  (x)  and  associated 
eigenvector  xk. 

X|c  <  Xfc+i  (eigenvalues) 

*k  *k  3  1 
so  that 

0  1 

1  -  -  J51  (20) 

xi  (1) 

which  now  gives  the  AR  coefficients. 

4.  The  moving  average  component  to  the  ARMA  model  Is  determined  using  the 
AR  coefficients. 


P  P 

rs  (n)  =  l  I  ak  am  rx  (n+m-k) 
k=0  m=0 


(21) 


0<  n  <q 


with  rs  (n)  =  0  n  >  q 

so  that  the  moving  average  component  of  the  spectrum  is  given  by 


Sw  (ej“)  =  l  rs  (n)  e~^n 
n=-g 

To  find  the  zero's  zk  of  (ej“) 


(22) 


SMA  (eJu)  “  !b0!2  n  (l-zke-J“)  (1  -  7keJ“)  (23) 

where  we  require  |zk  |  <  1;  all  the  zero's  of  S^(eJt“)  will  be  complex  conjugate 
pairs.  With 


q 


l 

k=0 


•  l  q 

bke*J“k  =  bQ  n 
k=0 


(1 


zt,e 


■ju)) 


(24) 
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we  can  compare  coefficients  of  e"^*  to  determine  b^.  Thus  for  the 


rational  ARMA  model 


Sx  (ej“  )  » 


l  rs  (n)  e-J“  n 
n=-q 

1  +  a,  +....+  an  e“JP“  I2 


where  Sx  (e^“)  Is  the  spectrum. 

2.  For  spectral  models  using  only  a  finite  number  of  observation  so  that 
we  have  only  estimates  of  the  autocorrelation  function,  Cadzow  [4]  has  identified 
the  following  procedures. 

a.  Moving  average  model  -  We  again  use  Blackman-Tukey  procedure,  l.e., 

q 

S  (eJ“)  -  l  w(n)r  (n)e"J“n  (26) 


where 


1  N 

rx(n)  *  -  l  x  (k+n)  x  (k) 
N  k=l 


-q<  n  <  q 


rx(n)  =  0  outside 

2n/(N-l )  0  <  n  <  (N-l)/2 

w(n)=  2-2n/(N-l)  (N-l)/2  <  n  <  (N-l)  (27) 

0  otherwise 

w(n)  Is  known  as  the  Bartlett  window. 

Any  "suitable"  window  w(n)  can  be  used  to  weigh  rx(n);  for  the  present,  we  use 
w(n)=l. 

b.  Autoregressive  model  (AR)  - 

1.  Compute  R(1,j)  =  rx  (1-j)  for  1  <  1,  j  <  p  +  1 


1  N 

R(1,j)  =  - -  T  x  (k+1-j)x(k) 

N- 1 1-j I  k-1 


(28) 


2.  Ra  =  | b0 | 2ej 

where  bo  Is  selected  to  normalize  first  component  of  a,  i.e.,  ^(l)  =  1. 

3*  %  <eJ") 

c.  Autoregressive  moving  average  (ARM*)  model  -  For  an  ARMA  model  the  AR 
components  are  calculated  first,  as  follows: 

1.  Compute  the  autocorrelation  matrix  estimate  Ri  whose  elements  Ri  (1,j) 
are  given  by 

Rl  (1,j)=  rx  (q+l+1-j)  (30) 

for  1  <  1  <  t  and  l<j<p+lso  that  Ri  Is  a  tX(p+l)  matrix.  The  rank  of 
Rl  is  equal  to  the  min  (t,  p+1).  If  the  number  of  eigenvalues  are  computed 
for  Ri*  Ri  it  will  be  equal  to  min  (t,  p+1)  so  that  if  t>>  p+1  we  can  be 
assured  that  rank  RiRi=p+l.  Hence  p,  the  order  of  AR  component  of  our 
model,  is  determined.  The  value  of  t  will  be  bounded,  typically  selected,  to 
bep£t_<N  -  q  -  1. 

2.  Solve  for  _a  In  equation  R*  WR^  =  o  £i  where  a  is  a  normalizing 
constant  so  that  1st  component  of  aP  is  1  and  W  Is  a  diagonal  weighting  matrix 
which  is  taken  as  I  (identity)  where  Ri  is  unbiased. 

.  t  7 

3.  Ri  W  Ri ( 1 , j )  =  l  w(m)  ?  (q+m+1)  ?  (q+m+i-j)  (31) 

m=l 

for  1  <  1,  j  <  p  +  1  with  w(m)  correspopnding  to  diagonal  elements  of  diagonal 
weighting  matrix  W. 


4.  Compute  (R*  W  R^)"*  and 
solve  for 


(Ri  “  Ri)_1H 


i°(i) 


which  normalizes  ji  so  that  1st  component  of  a0  is  1. 


|A  (eJ  “)|2  =  |  1  +  I  aP  e-J  J<|2 
p  ksl 


(32) 


(33) 


Moving  average  model.  Several  methods  exist  for  computing  the  moving  average 
coefficients  which  require  spectral  factorization  of  the  moving  average  spectrum. 

Method  #1  -  The  system  to  be  described  generates  a  moving  average  "residual 
signal"  obtained  by  passing  the  observed  data  x(t)  through  a  filter  whose 
transfer  function  Is  the  denominator  of  the  ARMA  model,  hence  an  AR(p),  thus 
producing  a  moving  average  output  whose  MA  spectrum  Bq  corresponds  to  the 
original  ARMA  model.  System  is  described  in  Figure  2. 

Residual  signal 

- - — * 

$ 


Figure  2. 

p+1  <  n  <  N  (34) 


P  7 

Sb(n)  ■  l  aic  x(n+k) 
k=0 


P  . 

Define  sf(n)  *  £  a^  x(n-k) 
k=0 


White  noise 

e  ■  ) 


Bq 

Observed  Data 

Ap 

X 

Ap 

1<  n  <  N-p 


rs(n)  = 


(N-p-n)  k*l 


l  sf(n+p+k)sf(p+k)+st>(n+k)sb(k)  (35) 


for  0<n<q  and  r$  (-n)  *  rs(n).  It  should  be  noted  that  for  a  moving  average 


model  rs(n)aQ  for  n>q+l.  The  spectrum  (ej“)  is  given  as  follows 


S^(eJ“)  =  |Bq  (e^MI2  *  l  w(n)  r$(n)e" 

n=-q 


jo>n 


(36) 


where  w(n)  = i 


N-p-n  \  / q+1  -  |n| 


(37) 


N-P 


q  +  1 


For  w=l 


|Bq(e^)|2  =  a2  [l+rs(l)  (z-4z)+. . .+r5( (38) 


where  z  =  eJ« 

In  general  we  can  find  the  zero's  zk  of  lBq(eJw)l2  *  Bq  (e^*0 )  !Tq(eJ“) 
and  carry  out  a  spectral  factorization 


|B  (e>)|2  «  b§  n  (l-zke->)  (l-zk  e>) 
H  k=l 


(39) 


so  that 


Bq(e>) 


b0  n  (l-zk  e’j“) 
k-1 


(40) 


Since 


Bq(e>)  -  l  bk  e 
H  k=0 


(41) 


we  have 


J  bk  e-j“k  ■  b0  u  (l-zke-J“) 


(42) 


By  equating  coefficients  we  determine  the  values  of  bk. 

Method  #2.  Assuming  the  AR  parameters  a^  are  known,  we  define 
function  rx( n)  and  Its  causal  image  rx(n)  as  follows 

rx(n)  *  E  {  x  (n+m)  x*(m)  }  i  n“°»  l1*  ±2»  ••• 

1 

rx(n)  =  rx(n)u(n)  -  -  rx  (0)  6  (n) 


where 


u(n)  *  standard  unit  step  * 


1  n>0 

0  n<0 


6(n)  ■  Kronecker  delta  sequence 

1  n  =  0 

0  n  *  0 


and  rx(-n)  =  rx  ^ 

Noting  that  rx(n)  =  rx(n)  +  rx  (-n)* 

we  have  Sx(eJ“)  =  Sx  {eJ“)  +  Sx  (e^) 

=  2  Re  [Sj  (eJ“)] 

which  Is  the  spectrum  of  the  signal 

P 

and  c(n)  =  rj(n)  +  £  ak  rj(n-k)  0  < 

=0  n  outside  interval 

Cs(e^“)  =  l  c(n)  e-jwn 
n=0 


Define 


where  s  «  max  (q,p) 


so  that 


Cs(e*°) 


=  [1+  l  an  e-J“n3  S+  (eJ“) 


Ap(e^“)  A*(< 

;j(0) 

A*p(ej“)  (Cs(e^)  +  ( 

:*{e^“)  Ap(ej“) 

Ap  (eJ“)  A*  (eJ“) 


However  Bqte^jBqteJ")  =  Aj(e**>)Cs(ej“)  +  Ap(e^)C*  (e*-) 

so  that  If  we  find  the  zero's  of  Sx  {ej“)  identified  as  zk  and  using  spectral 

factorization  we  can  determine  the  bk  by  equating  coefficients  in 


Q  q 

I  bk  e">k  =  b0^n^  (1-z^-1) 


for  |zkl<  1  (48) 


to  maintain  minimum  phase. 

.  ORDER  DETERMINATION 

A  fundamental  issue  in  applying  the  methods  that  have  been  presented  is 
that  the  order  of  the  model  needs  to  be  determined.  Methods  that  have  been 
used  include: 

1.  Levinson-Durbin  Method  for  pure  AR. 

2.  Test  for  autocorrelation  function,  rs(n)=0  with  n>q+l  for  pure  MA. 

A  new  method  receiving  considerable  study  is  the  singular  value  decomposition 
(SVD)  using  the  Frobenious  norm  of  a  mxn  matrix  difference,  A-B, 
defined  as 


1 1 A  -  B  II 


la,-i  -  b,-4  T 


The  SVD  Is  carried  out  as  follows: 


1.  Set  pe  »  p,  qe  »  q  for  t  »  p 

2.  Form  Re: 


Re  - 


rx(qe+l) 


rx(qe+t-D 


rx(qe"Pe+1J 

rx(qe-pe+2) 


rx(qe-pe+t) 


This  is  a  tx(pe+l)  matrix  which  satisfies 

(50) 

for  which  Re  =  U  E  V* 

where  U  and  V  are  unitary  matrices 

3.  Determine  eigenvalues  Akfe  of  RgR* 

4.  Take  akk  =+  /  x kk ,  called  the  singular  values,  which  are  ordered 
oll>a22>***>  ohh>0  where  h  =  min  (t,  pe+l) 

5.  Form  the  matrix  ReRe*  which  is  nonnegative  hermitian.  Using  the 
Gram- Schmidt  method  for  calculating  eigenvectors,  determine  the  columns  of  U, 
txt  matrix,  corresponding  to  ordered  orthonormal  eigenvectors  of  RgRg. 

6.  The  columns  of  V  are  the  orthonormal  eigenvectors  of  RgRg.  V  is  a 
(pe+l)x(pe+l)  matrix. 

7.  Etx(pe+1)  a  matrix  whose  elements  are  zero  except  possibly  along 
main  diagonal,  akk- 

Form  A(k)  =  U  lk  V* 

where  I  k  Is  obtained  by  setting  to  zero  all  but  its  "k"  largest  singular  val 


8. 


Finally,  we  define 


v(k)  ■ 


MA^ll 

!|A|| 


[gl \*  °2$+-"+  qk^1/2 

[®1?+  <}£+...+  ahh2]1/2 


(51) 


so  that  as  v(k) — >  1  as  k  — *  p,  the  order  of  our  model. 

The  above  approach  requires  that  singular  values  be  generated  for  carrying  out 
SVD.  Improved  accuracy  for  computing  singular  values  has  been  developed  by 
Dongarra  [8]. 

The  work  of  S.  Y.  Kung  [9]  using  state  space  variables  with  SVD  has  recently 
demonstrated  superior  spectral  estimates  than  estimates  obtained  in  terms  of 
transfer  function  parameters.  The  work  in  this  area  should  be  explored  for 
improved  estimates  for  perturbed  covariance  data. 

5.  ADAPTIVE  MODELING 

The  goal  of  this  investigation  is  to  develop  refined  models  which  adapt  to 
a  changing  environment.  Naturally,  it  becomes  appropriate  to  develop  adaptive 
models  which  conti nously  update  the  parameters  as  new  time  series  observations 
become  available  (x(N+l),  x  (N+2),...).  To  start  this  investigation  it  will  be 
appropriate  to  develop  a  class  of  adaptive  autocorrelations  estimates  and 
adaptive  algorithms  which  from  experience  have  proved  most  successful.  We  will 
define  an  adaptive  class  of  autocorrelation  estimators  as  follows. 


R(i.j) 


1 

(N  +  k2  -  ki) 


n+k2-l  _ 

l  x  ( k+l-i ) x( k+1- j ) 
k=ki 


(52) 


where  N  is  the  total  set  of  observations 


1  <  1  <  p  +  1 

1  <  j  <  p  +  1 
l  <  ki,  k2  <  p  +1 

with  Iq  and  k2  selected  so  that  the  number  of  lag  products  (N  +  k2  -  k^) 
is  at  least  equal  to  p+1. 

The  data  matrix  R  which  uses  R  (i,  j)  as  its  (i,  j)  elements  can  be 
written  in  the  form 


N  +  k2-ki 


xn 


where  XN  is  a  (N  +  k2  -ki)  x  (p  +  1)  data  matrix  such  that 
XN  (i,  j)  ■  x  (ki  +  i  -  j) 
and 

1  <  i<  N  +  k2  -  ki 
1  <  j  <  p  +  1 

It  should  be  noted  that  we  set  x  =  0  whenever  kj  +i-j  falls  outside  the 
observation  set  1  <  n  <  N. 

We  have  now  established  the  data  matrix  to  be  used  in  adaptive  modeling 

A 

Experience  has  shown  that  R  is  unbiased  and  provides  consistent  statistics 
when  ki=p+l  and  k2al.  This  method  is  called  the  covariance  method. 

AR  Model  -  Covariance  Method: 

As  indicated,  with  k^  *  p+1  and  k2=l» 

*N  XN  *N  *  (N+k2-ki)  |b0l2  ej  (54) 

where  bo  is  selected  to  normalize  1st  component  of  to  1. 

Now  let 


XN+1  XN+1  s  XN  XN  +  iN+l  £n+1  p+1^ 


so  that 


where  xjtj+i  =  [x  (N+l),  x  { N ) ,  ...,x  (N+l-p)] 

For  k2=l  and  N  =  p+1,  and  the  initial  value  xnxn=xp+i  Xp+^ 


*  p+1  _ 

Vl  xpH  3  S  x  (kn'1)  x  {k+1-J)  (56) 

k=ki 


for  l<i<p+l  and  l<j<p+l.  As  can  be  seen  from  the  expression  for  ajj+i  we  need 
to  compute  [XN+1  Xn+i3"a .  This  is  given  by 

(2n+1  ^N+l  ^ 


Bh+iW1  *  - 


d  +iN+l  *N*+1> 


(57) 


where 


&I+1  =  *  N+l  for  N>P+kl 

Accordingly,  using  Gaussian  elimination,  we  can  calculate  [X^X^”1  for  N^p+kj 

"it  i 

so  that  for  all  N>  p+klt  we  can  calculate  CXN+1XN+i] “A  using  the  above  expression. 
This  will  be  used  for  updating  the  parameter  ajj+i . 

AR  adaptive  algorithms  then  require  the  following  steps: 

Step  1:  Input  data:  x(N+l),  [X^XN]-1 
Step  2:  Compute:  Exn+ixn+i]_1 
Step  3:  Let  £  »  ^n+^N+I^"1 

Step  4:  a^  =  c(l)"*  £  where  c(l)  is  the  1st  component  of  £ 

The  problem  of  AR  order  determination  still  remains  as  in  the  non-adaptive  case. 
The  approach  recommended  In  using  raw  time  series  data  where  R  (i,j)  =rx  ( i - j ) 
for  l<i,j<p+l  is  to  find  the  order  "pi"  for  which  R  has  ( p-pi )  of  its  eigenvalues 
sufficiently  close  to  zero  for  all  p>pi»  This  could  be  carried  out  using  SVD. 

ARMA  Adaptive  Modeling: 

A 

As  in  the  AR  adaptive  model  we  need  to  define  Rj  (i,  j)  which  we  do  by  using 


Rl(i  J) 


1  N-q+2+k2_ 

■  l  x  (k+l-i)x(k+q+2-j ) 

(N+k2-ki-q-l)  k=ki 


(58) 


with  1<1 <t  and  l<j<p+l . 


Rl(1,j)  a  rx  (q+l+1-j) 


The  number  of  lag  products,  (N+k2-ki~q-l) ,  Is  selected  such  that 
(N*k2-ki-q-l)>p+l  and  1  <  ki  <  t  1  <  k2  <  P+1. 

As  In  the  AR  adaptive  model  the  covariance  method  had  preferred  properties.  In 
the  ARMA  adaptive  model  the  covariance  method  with  kj=t  amd  k2=l  has  similarly 
been  demonstrated  as  being  unbiased  and  statistically  consistent. 

ARMA  Covariance  Method: 

As  Indicated  we  constrain  the  parameters  k^,  k2  as  shown:  kj  =  t,  k2  =  1, 
then  we  write  Ri  as  follows 


where 

for 


1 


N+k2-ki-q-l 


1  Y3  xN 


XM(1,j)  *  x  (k!+q+l+1-j) 

1  <  1  <  N+k2  -  ki  -  q  -  1 
1  <  j  <  p+1 


(59) 


wl  th 

and 


with 


YN  (1,j)  =  x  (ki+1-j) 

1  <  1  <  N+k2  -  kj  -  q  -  1 
1  <  j  <  t 

x(n)  =  0  for  n>N  or  n<l 


To  determine  the  AR  parameters  we  need  to  solve 

XN  YN  YN  XN  iN  b  ail  (60) 

where  a  is  selected  for  the  first  component  of  a#  equal  to  1. 

When  the  updated  a^+j  parameter  Is  calculated  we  will  require  • 


This  can  be  calculated  using 


yN+1xN+1  =  YNXN  +  2n2n  {N>t) 

where  XN  *  Cx(N+l),  x(N),...,x{N+l-p)] 

=  Cx(M-q) ,x(N-q-l) .... ,x(N+l-q-t)] 
so  that  for  N>t 


xN+1y 


N+1yN+1] 


‘N+l 


-  xnynynxn 


hi  hi 


h\  In 


hi  hi 


(61) 


and 

hi  =  &  yNxN 

The  best  overall  ARMA  adaptive  model  performance  is  the  covariance  method 
for  k2=l  and  ki=t  (ki  can  range  In  the  Interval  1  <  k i  <  t);  the  ARMA  adaptive 
model  algorithm  for  k2“l  and  l<ki<t1s  given  as  follows: 

Step  0:  The  Input  to  commence  the  algorithm  atN  =  q+  p  +  ki  +  l 
is  xN*  xN  and  [xj  rM  rj  xNr‘ 

which  can  be  calculated  by  Gaussian  elimination. 

Step  1:  N  =  q+  p  +  ki  +  l 
Step  2:  Compute  YjJj+j  XN+j  from 

yN+1  xN+1  =  yN  xN  +  &l  JSN 

xjj  =  [x(N+l )  ,x(N) , . . .  ,x(N+l-p)] 

=  Cx(N-q) ,x(N-q-l ) ,. . . ,x( N+l-q-t) ] 

Step  3:  in  =  Zulu  h 
Step  4:  =  zjj 

n  =  b 

‘I1  -  (xn’,»vnxn)'1 

Compute  [Aj  +  u£  vj]"*  from 


28 


[A  +  _u*v]“l  =  A“1 


[A-lu*]  [vA-l] 
(l+vA"l  jj*) 


Step  5:  ug  *  -£n»  *  In 

Ag1  »  [Aj  +  Ui  vj]-1 

Compute  [A£  +  U2  ^3"^  from  step  4 
Step  6:  U3  =  (^J) 

n  3  IN 

A31  =  [A2  +  uj  V2]"1 

Compute  [A3  +  .U3  V3]"1  *  CxJ+^N+^N+^N+l^*1 


from  step  4. 

Step  7: 

1  *  [xn+iyn+iyn+ixn+i3-1  ll 

where 

afl+1  ■  [c(l)]-1  c 

c(l)  is  the  first  component  of  £. 

Step  8: 

Let  N  =  N  +  1,  go  to  step  2. 

We  have  presented  a  sufficient  number  of  useful  methods  that  should  be 
explored  in  stochastic  modeling  for  any  application  you  may  have  in  mind.  The 
application  of  SVO  provides  an  important  method  for  determining  order  of  the 
model.  The  recent  results  by  S.  Kung  using  state  space  and  SVD  appear  to 
provide  higher  resolution  models  and  should  be  evaluated  for  Army  applications. 
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